#ifndef _MESH_
#define _MESH_

#include "Domain.h"
#include "Dofs.h"
#include "Grid.h"
#include <map>
#include <functional>
#include <fstream>
#include <sstream>
#include <set>
using namespace std;
/**
 * @brief 定义域划分成网格Mesh,Mesh包含小网格Grid,每个Grid的类型由TemplateElement确定
 *
 * @tparam int DIM 维度
 */

vector<int> GridId1toId2(int i, int _segx);
int DofsId2toId1(int i, int j, int _segx, int _order);

template <unsigned int DIM>
class Mesh
{
public:
  /**
   * @brief 默认构造函数
   */
  Mesh() = default;

  /**
   * @brief 默认析构函数
   */
  virtual ~Mesh() = default;
  /**
   * @brief 返回自由度的总个数
   * @details _GridList.size()*getNumDofs();
   * @return
   */
  virtual int getTotalNumDofs() = 0;
  /**
   * @brief 返回网格的总个数
   * @details _GridList.size()
   * @return [description]
   */
  virtual int getNumGrid() = 0;
  /**
   * @brief 返回存储Grid的vector
   * @details _GridList
   * @return [description]
   */
  virtual vector<Grid<DIM>> &getGrid() = 0;
  /**
   * @brief 返回序号为_index的Grid
   *
   * @param _index Grid的序号
   * @return _GridList[i];
   */
  virtual Grid<DIM> &getGrid(int _index) = 0;
  /**
   * @brief   返回vector，存储边界上的自由度
   * @details 对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
   *
   * @param component 试图对边界分段，并标号
   * @return 存储边界上自由度的vector
   */
  virtual vector<Grid<2>> getBoundaryGrid() = 0;
  // virtual void setBoundaryIndex() = 0;
  virtual set<Dofs<2>> getBoundaryDofs() = 0;
  virtual vector<pair<int, int>> getIndexofGrid(int i) = 0;
  virtual vector<int> &getSegment() = 0;
};

class Q1Mesh : public Mesh<2>
{
public:
  Q1Mesh() = default;
  /**
   * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
   * @param domain 定义域
   * @param seg 每个维度要划分的段数
   */
  Q1Mesh(Domain<2> *_domain, vector<int> _seg);
  ~Q1Mesh() = default;

  int getTotalNumDofs();
  int getNumGrid();
  vector<Grid<2>> &getGrid();
  Grid<2> &getGrid(int i);
  vector<Grid<2>> getBoundaryGrid();
  vector<int> &getSegment();
  // void setBoundaryIndex();
  set<Dofs<2>> getBoundaryDofs();
  vector<pair<int, int>> getIndexofGrid(int i);//取名取得不好

protected:
  vector<Grid<2>> gridList;           /*< 存放所有的Grid*/
  vector<pair<int, int>> indexonMesh; /*< 存放mesh上的index,dofsindex,boundarymark */
  vector<int> seg;                    /*< 存放各维度分段数*/
};

Q1Mesh::Q1Mesh(Domain<2> *_domain, vector<int> _seg)
{
  assert(_seg.size() == 2);
  Point<2> p0 = _domain->getVertex(0);
  Point<2> p2 = _domain->getVertex(2);
  double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
  double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
  for (int j = 0; j < _seg[1]; ++j)
  {
    for (int i = 0; i < _seg[0]; ++i)
    {
      vector<Point<2>> pnts;
      pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
      pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
      Grid<2> grid(pnts);
      gridList.push_back(grid);
    }
  }
  for (int i = 0; i < (_seg[0] + 1) * (_seg[1] + 1); i++)
  {
    indexonMesh.push_back({i, 0});
  }
  seg = _seg;
}

int Q1Mesh::getNumGrid()
{
  return gridList.size();
}

Grid<2> &Q1Mesh::getGrid(int i)
{
  return gridList[i];
}

vector<Grid<2>> &Q1Mesh::getGrid()
{
  return gridList;
}

vector<Grid<2>> Q1Mesh::getBoundaryGrid()
{
  vector<Grid<2>> boundary_grid;
  int k = 0;
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  k = seg[0];
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] - 1;
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] * (seg[1] - 1);
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  return boundary_grid;
}

vector<int> &Q1Mesh::getSegment()
{
  return seg;
}

vector<pair<int, int>> Q1Mesh::getIndexofGrid(int i)
{
  vector<pair<int, int>> idx;
  int a = seg[0];
  idx.push_back(indexonMesh[i + i / a]);
  idx.push_back(indexonMesh[i + i / a + 1]);
  idx.push_back(indexonMesh[i + i / a + a + 2]);
  idx.push_back(indexonMesh[i + i / a + a + 1]);
  return idx;
}

// void Q1Mesh::setBoundaryIndex()
// {
//     vector<int> k = {0, seg[0] + 1, seg[0], (seg[0] + 1) * seg[1] + 1};
//     for (int i = 0; i < seg[0]; i++)
//     {
//         indexonMesh[k[0]].second = 1;
//         k[0]++;
//         indexonMesh[k[3]].second = 1;
//         k[3]++;
//     }
//     for (int i = 0; i < seg[1]; i++)
//     {
//         indexonMesh[k[1]].second = 1;
//         k[1] += seg[0] + 1;
//         indexonMesh[k[2]].second = 1;
//         k[2] += seg[0] + 1;
//     }
// }

set<Dofs<2>> Q1Mesh::getBoundaryDofs()
{
  set<Dofs<2>> boundaryDofs;
  vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
  vector<int> m = {0, seg[0] + 1, seg[0], (seg[0] + 1) * seg[1] + 1};
  for (int i = 0; i < seg[0]; i++)
  {
    Dofs<2> d0 = {gridList[k[0]].getVertices(0), m[0]};
    Dofs<2> d3 = {gridList[k[3]].getVertices(2), m[3]};
    boundaryDofs.insert(d0);
    boundaryDofs.insert(d3);
    k[0]++;
    k[3]++;
    m[0]++;
    m[3]++;
  }
  for (int i = 0; i < seg[1]; i++)
  {
    Dofs<2> d1 = {gridList[k[1]].getVertices(3), m[1]};
    Dofs<2> d2 = {gridList[k[2]].getVertices(1), m[2]};
    boundaryDofs.insert(d1);
    boundaryDofs.insert(d2);
    k[1] += seg[0];
    k[2] += seg[0];
    m[1] += seg[0] + 1;
    m[2] += seg[0] + 1;
  }
  return boundaryDofs;
}

int Q1Mesh::getTotalNumDofs()
{
  return (seg[0] + 1) * (seg[1] + 1);
}

class Q2Mesh : public Mesh<2>
{
public:
  Q2Mesh() = default;
  /**
   * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
   * @param domain 定义域
   * @param seg 每个维度要划分的段数
   */
  Q2Mesh(Domain<2> *_domain, vector<int> _seg);
  ~Q2Mesh() = default;

  int getTotalNumDofs();
  int getNumGrid();
  vector<Grid<2>> &getGrid();
  Grid<2> &getGrid(int i);
  vector<Grid<2>> getBoundaryGrid();
  vector<int> &getSegment();
  // void setBoundaryIndex();
  set<Dofs<2>> getBoundaryDofs();
  vector<pair<int, int>> getIndexofGrid(int i);

protected:
  vector<Grid<2>> gridList;           /*< 存放所有的Grid*/
  vector<pair<int, int>> indexonMesh; /*< 存放mesh上的index */
  vector<int> seg;                    /*< 存放各维度分段数*/
};

Q2Mesh::Q2Mesh(Domain<2> *_domain, vector<int> _seg)
{
  assert(_seg.size() == 2);
  Point<2> p0 = _domain->getVertex(0);
  Point<2> p2 = _domain->getVertex(2);
  double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
  double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
  for (int j = 0; j < _seg[1]; ++j)
  {
    for (int i = 0; i < _seg[0]; ++i)
    {
      vector<Point<2>> pnts;
      pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
      pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
      Grid<2> grid(pnts);
      gridList.push_back(grid);
    }
  }
  for (int i = 0; i < (2 * _seg[0] + 1) * (2 * _seg[1] + 1); i++)
  {
    indexonMesh.push_back({i, 0});
  }
  seg = _seg;
}

int Q2Mesh::getNumGrid()
{
  return gridList.size();
}

Grid<2> &Q2Mesh::getGrid(int i)
{
  return gridList[i];
}

vector<Grid<2>> &Q2Mesh::getGrid()
{
  return gridList;
}

vector<int> &Q2Mesh::getSegment()
{
  return seg;
}

int Q2Mesh::getTotalNumDofs()
{
  return (2 * seg[0] + 1) * (2 * seg[1] + 1);
}

// void Q2Mesh::setBoundaryIndex()
// {
//     vector<int> k = {0, 2 * seg[0] + 1, 2 * seg[0], (2 * seg[0] + 1) * (2 * seg[1]) + 1};
//     for (int i = 0; i < 2 * seg[0]; i++)
//     {
//         indexonMesh[k[0]].second = 1;
//         k[0]++;
//         indexonMesh[k[3]].second = 1;
//         k[3]++;
//     }
//     for (int i = 0; i < 2 * seg[1]; i++)
//     {
//         indexonMesh[k[1]].second = 1;
//         k[1] += 2 * seg[0] + 1;
//         indexonMesh[k[2]].second = 1;
//         k[2] += 2 * seg[0] + 1;
//     }
// }

vector<pair<int, int>> Q2Mesh::getIndexofGrid(int i)
{
  vector<pair<int, int>> idx;
  int a = seg[0];
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2]);
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2 + 1]);
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + i % a * 2 + 2]);
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2 + 2]);
  idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2 + 2]);
  idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2 + 1]);
  idx.push_back(indexonMesh[(i / a + 1) * (4 * a + 2) + i % a * 2]);
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2]);
  idx.push_back(indexonMesh[i / a * (4 * a + 2) + 2 * a + 1 + i % a * 2 + 1]);
  return idx;
}

vector<Grid<2>> Q2Mesh::getBoundaryGrid()
{
  vector<Grid<2>> boundary_grid;
  int k = 0;
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  k = seg[0];
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] - 1;
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] * (seg[1] - 1);
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  return boundary_grid;
}

set<Dofs<2>> Q2Mesh::getBoundaryDofs()
{
  set<Dofs<2>> boundaryDofs;
  vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
  vector<int> m = {0, 2 * seg[0] + 1, 2 * seg[0], (2 * seg[0] + 1) * 2 * seg[1] + 1};
  for (int i = 0; i < seg[0]; i++)
  {
    Point<2> p0_0 = gridList[k[0]].getVertices(0);
    Point<2> p0_1 = gridList[k[0]].getVertices(1);
    Dofs<2> d0_0 = {p0_0, m[0]};
    Dofs<2> d0_1 = {(p0_0 + p0_1) / 2, m[0] + 1};
    Point<2> p2_2 = gridList[k[3]].getVertices(2);
    Point<2> p2_3 = gridList[k[3]].getVertices(3);
    Dofs<2> d2_5 = {(p2_2 + p2_3) / 2, m[3]};
    Dofs<2> d2_4 = {p2_2, m[3] + 1};
    boundaryDofs.insert(d0_0);
    boundaryDofs.insert(d0_1);
    boundaryDofs.insert(d2_5);
    boundaryDofs.insert(d2_4);
    k[0]++;
    k[3]++;
    m[0] += 2;
    m[3] += 2;
  }
  for (int i = 0; i < seg[1]; i++)
  {
    Point<2> p3_0 = gridList[k[1]].getVertices(0);
    Point<2> p3_3 = gridList[k[1]].getVertices(3);
    Dofs<2> d3_7 = {(p3_0 + p3_3) / 2, m[1]};
    Dofs<2> d3_6 = {p3_3, m[1] + 2 * seg[0] + 1};
    Point<2> p1_1 = gridList[k[2]].getVertices(1);
    Point<2> p1_2 = gridList[k[2]].getVertices(2);
    Dofs<2> d1_2 = {p1_1, m[2]};
    Dofs<2> d1_3 = {(p1_1 + p1_2) / 2, m[2] + 2 * seg[0] + 1};
    boundaryDofs.insert(d3_7);
    boundaryDofs.insert(d3_6);
    boundaryDofs.insert(d1_2);
    boundaryDofs.insert(d1_3);
    k[1] += seg[0];
    k[2] += seg[0];
    m[1] += 2 * (2 * seg[0] + 1);
    m[2] += 2 * (2 * seg[0] + 1);
  }
  return boundaryDofs;
}

class Q4Mesh : public Mesh<2>
{
public:
  Q4Mesh() = default;
  /**
   * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
   * @param domain 定义域
   * @param seg 每个维度要划分的段数
   */
  Q4Mesh(Domain<2> *_domain, vector<int> _seg);
  ~Q4Mesh() = default;

  int getTotalNumDofs();
  int getNumGrid();
  vector<Grid<2>> &getGrid();
  Grid<2> &getGrid(int i);
  vector<Grid<2>> getBoundaryGrid();
  vector<int> &getSegment();
  // void setBoundaryIndex();
  set<Dofs<2>> getBoundaryDofs();
  vector<pair<int, int>> getIndexofGrid(int i);

protected:
  vector<Grid<2>> gridList;           /*< 存放所有的Grid*/
  vector<pair<int, int>> indexonMesh; /*< 存放mesh上的index */
  vector<int> seg;                    /*< 存放各维度分段数*/
};

Q4Mesh::Q4Mesh(Domain<2> *_domain, vector<int> _seg)
{
  assert(_seg.size() == 2);
  Point<2> p0 = _domain->getVertex(0);
  Point<2> p2 = _domain->getVertex(2);
  double xLen = (p2[0] - p0[0]) / (1.0 * _seg[0]);
  double yLen = (p2[1] - p0[1]) / (1.0 * _seg[1]);
  for (int j = 0; j < _seg[1]; ++j)
  {
    for (int i = 0; i < _seg[0]; ++i)
    {
      vector<Point<2>> pnts;
      pnts.push_back({p0[0] + i * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + j * yLen});
      pnts.push_back({p0[0] + (i + 1) * xLen, p0[1] + (j + 1) * yLen});
      pnts.push_back({p0[0] + i * xLen, p0[1] + (j + 1) * yLen});
      Grid<2> grid(pnts);
      gridList.push_back(grid);
    }
  }
  for (int i = 0; i < (4 * _seg[0] + 1) * (4 * _seg[1] + 1); i++)
  {
    indexonMesh.push_back({i, 0});
  }
  seg = _seg;
}

int Q4Mesh::getNumGrid()
{
  return gridList.size();
}

Grid<2> &Q4Mesh::getGrid(int i)
{
  return gridList[i];
}

vector<Grid<2>> &Q4Mesh::getGrid()
{
  return gridList;
}

vector<int> &Q4Mesh::getSegment()
{
  return seg;
}

int Q4Mesh::getTotalNumDofs()
{
  return (4 * seg[0] + 1) * (4 * seg[1] + 1);
}

vector<pair<int, int>> Q4Mesh::getIndexofGrid(int i)
{
  vector<pair<int, int>> idx;
  int x = GridId1toId2(i, seg[0])[0];
  int y = GridId1toId2(i, seg[0])[1];
  for (int j = 0; j < 5; j++)
    for (int i = 0; i < 5; i++)
      idx.push_back(indexonMesh[DofsId2toId1(4 * x + i, 4 * y + j, seg[0], 4)]);
  return idx;
}

vector<Grid<2>> Q4Mesh::getBoundaryGrid()
{
  vector<Grid<2>> boundary_grid;
  int k = 0;
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  k = seg[0];
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] - 1;
  for (int i = 0; i < seg[1] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k += seg[0];
  }
  k = seg[0] * (seg[1] - 1);
  for (int i = 0; i < seg[0] - 1; i++)
  {
    boundary_grid.push_back(gridList[k]);
    k++;
  }
  return boundary_grid;
}

set<Dofs<2>> Q4Mesh::getBoundaryDofs()
{
  set<Dofs<2>> boundaryDofs;
  vector<int> k = {0, 0, seg[0] - 1, seg[0] * (seg[1] - 1)};
  vector<int> m = {0, 4 * seg[0] + 1, 4 * seg[0], (4 * seg[0] + 1) * 4 * seg[1] + 1};
  for (int i = 0; i < seg[0]; i++)
  {
    Point<2> p0_0 = gridList[k[0]].getVertices(0);
    Point<2> p0_1 = gridList[k[0]].getVertices(1);
    Dofs<2> d0_0 = {p0_0, m[0]};
    Dofs<2> d0_1 = {p0_0 * 0.75 + p0_1 * 0.25, m[0] + 1};
    Dofs<2> d0_2 = {p0_0 * 0.5 + p0_1 * 0.5, m[0] + 2};
    Dofs<2> d0_3 = {p0_0 * 0.25 + p0_1 * 0.75, m[0] + 3};
    Point<2> p2_2 = gridList[k[3]].getVertices(2);
    Point<2> p2_3 = gridList[k[3]].getVertices(3);
    Dofs<2> d2_21 = {p2_3 * 0.75 + p2_2 * 0.25, m[3]};
    Dofs<2> d2_22 = {p2_3 * 0.5 + p2_2 * 0.5, m[3] + 1};
    Dofs<2> d2_23 = {p2_3 * 0.25 + p2_2 * 0.75, m[3] + 2};
    Dofs<2> d2_24 = {p2_2, m[3] + 3};
    boundaryDofs.insert(d0_0);
    boundaryDofs.insert(d0_1);
    boundaryDofs.insert(d0_2);
    boundaryDofs.insert(d0_3);
    boundaryDofs.insert(d2_21);
    boundaryDofs.insert(d2_22);
    boundaryDofs.insert(d2_23);
    boundaryDofs.insert(d2_24);
    k[0]++;
    k[3]++;
    m[0] += 4;
    m[3] += 4;
  }
  for (int i = 0; i < seg[1]; i++)
  {
    Point<2> p3_0 = gridList[k[1]].getVertices(0);
    Point<2> p3_3 = gridList[k[1]].getVertices(3);
    Dofs<2> d3_5 = {p3_0 * 0.75 + p3_3 * 0.25, m[1]};
    Dofs<2> d3_10 = {p3_0 * 0.5 + p3_3 * 0.5, m[1] + 4 * seg[0] + 1};
    Dofs<2> d3_15 = {p3_0 * 0.25 + p3_3 * 0.75, m[1] + 2 * (4 * seg[0] + 1)};
    Dofs<2> d3_20 = {p3_3, m[1] + 3 * (4 * seg[0] + 1)};
    Point<2> p1_1 = gridList[k[2]].getVertices(1);
    Point<2> p1_2 = gridList[k[2]].getVertices(2);
    Dofs<2> d1_4 = {p1_1, m[2]};
    Dofs<2> d1_9 = {p1_1 * 0.75 + p1_2 * 0.25, m[2] + 4 * seg[0] + 1};
    Dofs<2> d1_14 = {p1_1 * 0.5 + p1_2 * 0.5, m[2] + 2 * (4 * seg[0] + 1)};
    Dofs<2> d1_19 = {p1_1 * 0.25 + p1_2 * 0.75, m[2] + 3 * (4 * seg[0] + 1)};
    boundaryDofs.insert(d3_5);
    boundaryDofs.insert(d3_10);
    boundaryDofs.insert(d3_15);
    boundaryDofs.insert(d3_20);
    boundaryDofs.insert(d1_4);
    boundaryDofs.insert(d1_9);
    boundaryDofs.insert(d1_14);
    boundaryDofs.insert(d1_19);
    k[1] += seg[0];
    k[2] += seg[0];
    m[1] += 4 * (4 * seg[0] + 1);
    m[2] += 4 * (4 * seg[0] + 1);
  }
  return boundaryDofs;
}

vector<int> GridId1toId2(int i, int _segx)
{
  vector<int> a = {i % _segx, i / _segx};
  return a;
}

int DofsId2toId1(int i, int j, int _segx, int _order)
{
  return i + (_order * _segx + 1) * j;
}

// template<unsigned int DIM>
// class UnstructedMesh: public Mesh<DIM>
// {
// public:
//     // virtual void readMesh();
//     // virtual void writeMesh();
//     // virtual void distributeDofs(int _order_of_element);
//     virtual void InputDofData() = 0;
//     virtual void InputEdgeData() = 0;
//     virtual void InputEleData() = 0;
// };

// class EasyMesh :public Mesh<2>
// {
//   private:
//     vector<pair<int, int>> indexonMesh; /*< 存放mesh上的dofs:index,boundarymark */
//     string meshName;
//     long int nDofs;
//     long int nGrids;
//     long int nEdges;
//     long int nBoundaryEdges;
//     vector<Dofs<2> > dofsList;
//     vector<Dofs<2> > boundaryDofs;
//     vector<Grid<2> > gridsList;           /*< 存放所有的Grid*/

//     vector<vector<int>> boundaryEdges;
//     // not use function
//     vector<int> seg={-1,-1};
//     virtual vector<int>& getSegment(){return seg; };

//   public:
//     EasyMesh();
//     EasyMesh(string meshname);
//     void inputDofData();
//     void inputEdgeData();
//     void inputGridData();
//     virtual int getNumGrid(){return nGrids;};
//     virtual int getTotalNumDofs(){return nDofs;};
//     virtual vector<Grid<2> >& getGrid();
//     virtual vector<Grid<2> > getBoundaryGrid();
//     virtual set<Dofs<2> > getBoundaryDofs();
//     virtual vector<pair<int,int> > getIndexofGrid(int i);

//     //int getBndedgeNum(){return _nbndedge;};
//     //Dofs<2> getDofs(int i){return _DofsList[i];};
//     //vector<int> getElement(int i){return _ElementData[i];};
//     //vector<Dofs<2>>& getBoundary(){return _BndDofs;};
//     //bool IsBoundary(int index){return !(std::find(_BndDofs.begin(),_BndDofs.end(),index) == _BndDofs.end());};
//     //vector<int> getBndEdge(int i){return _BndEdges[i];};
// };

// EasyMesh::EasyMesh(string meshname)
// {
//   meshName = meshname;
//   inputDofData();
//   inputEdgeData();
//   inputGridData();
// }

// void EasyMesh::inputDofData()
// {
//   ifstream fin;
//   fin.open(meshName+".n");
//   string line;
//   stringstream ss;
//   getline(fin, line);
//   ss.clear();
//   ss.str(line);
//   ss >> nDofs >> nGrids >> nEdges;
//   dofsList.resize(nDofs);
//   gridsList.resize(nGrids);
//   int a, bndflag;
//   double x, y;
//   for (int i = 0; i < nDofs; i++)
//   {
//     getline(fin, line);
//     ss.clear();
//     ss.str(line);
//     ss >> a >> x >> y >> bndflag;

//     Point<2> temp_point({x, y});
//     Dofs<2> temp_dof(temp_point, i);
//     dofsList[i] = temp_dof;
//     if (bndflag)
//       boundaryDofs.push_back(temp_dof);
//   }
// }

// void EasyMesh::inputEdgeData()
// {
//   ifstream fin;
//   fin.open(meshName+".s");
//   string line;
//   stringstream ss;
//   getline(fin, line);
//   ss.clear();
//   ss.str(line);
//   //int _nedge;
//   //ss >> _nedge;
//   int a, d0, d1, bndflag;
//   for (int i = 0; i < nEdges; i++)
//   {
//     getline(fin, line);
//     ss.clear();
//     ss.str(line);
//     ss >> a >> d0 >> d1 >> a >> a >> bndflag;
//     if (bndflag)
//     {
//       vector<int> temp_edge(2);
//       temp_edge[0] = d0;
//       temp_edge[1] = d1;
//       boundaryEdges.push_back(temp_edge);
//     }
//   }
//   nBoundaryEdges = boundaryEdges.size();
// }

// void EasyMesh::inputGridData()
// {
//   ifstream fin;
//   fin.open(meshName+".e");
//   string line;
//   stringstream ss;
//   getline(fin, line);
//   //ss.clear();
//   //ss.str(line);
//   //ss >> nGrids;
//   //gridsList.resize(nGrids);

//   int a, d0, d1 ,d2;
//   for (int i = 0; i < nGrids; i++)
//   {
//     getline(fin, line);
//     ss.clear();
//     ss.str(line);
//     ss >> a >> d0 >> d1 >> d2;
//     vector<Point<2> > points= {dofsList[d0].getPoint(),dofsList[d1].getPoint(),dofsList[d2].getPoint()};
//     Grid<2> grid(points);
//     gridsList[i] = grid;
//     //gridList[i].resize(3);
//     //_GridData[i][0] = d0;
//     //_GridData[i][1] = d1;
//     //_GridData[i][2] = d2;
//   }
// }

// vector<Grid<2> >& EasyMesh::getGrid()
// {
//   set<Dofs<2> > s;
// }
#else
#endif
